Mendelian Randomization analysis of DeepSEA generated data

library(MendelianRandomization)

Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

data_dir = "../../dat"
seq_predictions_all = read.csv(file.path(data_dir, "means_and_uncertainties.csv"))
seq_predictions = subset(seq_predictions_all, seq_num == 1)

str(seq_predictions)

'data.frame':   3001 obs. of  6 variables:
 $ seq_num    : int  1 1 1 1 1 1 1 1 1 1 ...
 $ mut        : Factor w/ 3000 levels "0A","0C","0G",..: 2729 1 2 3 342 343 344 674 675 676 ...
 $ X_pred_mean: num  2.33e-03 3.50e-05 -1.89e-05 -1.57e-05 -1.33e-05 ...
 $ X_pred_var : num  0.0586 0.06 0.06 0.06 0.06 ...
 $ Y_pred_mean: num  1.71e-04 8.74e-06 -2.57e-06 -1.09e-06 -2.20e-06 ...
 $ Y_pred_var : num  0.00651 0.00777 0.00777 0.00777 0.00777 ...

bx <- unlist(seq_predictions["X_pred_mean"])
bxse <- unlist(seq_predictions["X_pred_var"])
by <- unlist(seq_predictions["Y_pred_mean"])
byse <- unlist(seq_predictions["Y_pred_var"])

MRInputObject <- mr_input(bx = bx, 
                          bxse = bxse, 
                          by = by, 
                          byse = byse)

cis = c()
for (seq in (1:50)) {
  seq_predictions = subset(seq_predictions_all, seq_num == seq)
  
  bx <- unlist(seq_predictions["X_pred_mean"])
  bxse <- unlist(seq_predictions["X_pred_var"])
  by <- unlist(seq_predictions["Y_pred_mean"])
  byse <- unlist(seq_predictions["Y_pred_var"])
  MRInputObject <- mr_input(bx = bx, 
                          bxse = bxse, 
                          by = by, 
                          byse = byse)
  EggerObject <- mr_egger(MRInputObject, robust = FALSE, alpha = .05)
  # cis = append(cis, c(EggerObject$CILower.Est, EggerObject$CIUpper.Est))
  print(c(EggerObject$CILower.Est, EggerObject$CIUpper.Est))
}

[1] 0.0590192 0.1159848
[1] 0.3233909 0.5597846
[1] 0.03235503 0.05877458
[1] 0.1188956 0.2318077
[1] 0.1645073 0.2648561
[1] 0.2289623 0.4052021
[1] 0.3443012 0.6641439
[1] 0.1756041 0.3033061
[1] 0.1583441 0.2399981
[1] 0.1889218 0.3005885
[1] 0.1827082 0.3066469
[1] 0.05220543 0.10226116
[1] 0.91425 1.57727
[1] 0.2863903 0.4708377
[1] 0.1051252 0.1822457
[1] 0.03974643 0.07074950
[1] 0.07505532 0.11982899
[1] 0.3706091 0.5714237
[1] 0.4507038 0.7233038
[1] 0.2518358 0.4329524
[1] 0.02264022 0.03560069
[1] 0.2323919 0.3507465
[1] 0.07638818 0.12772401
[1] 0.3171318 0.4806720
[1] 0.2739240 0.4173557
[1] 0.03117739 0.06000044
[1] 0.04368776 0.12557530
[1] 0.07505245 0.14740756
[1] 0.04114612 0.11159205
[1] 0.1269643 0.2866314
[1] 0.2601743 0.6519868
[1] 0.06622755 0.11109060
[1] 1.109005 1.907397
[1] 0.04254034 0.09048482
[1] 0.1153641 0.1847817
[1] 0.4348032 0.9299903
[1] 0.1471299 0.3256183
[1] 0.08097467 0.16127508
[1] 0.2559104 0.4957051
[1] 0.4693041 0.8698878
[1] 0.02719744 0.05512086
[1] 0.1165405 0.2815027
[1] 0.1989373 0.4144706
[1] 0.06628877 0.14424319
[1] 0.1635376 0.4410249
[1] 0.1459457 0.3396873
[1] 0.03371326 0.09627066
[1] 0.1937338 0.4058515
[1] 0.9633243 1.6609768
[1] 0.1370262 0.3841660

# dim(cis) <- c(2, 50)
# print(cis)

mr_plot(MRInputObject)

install.packages("devtools")

Installing package into ‘/home/stephenmalina/R/x86_64-pc-linux-gnu-library/3.6’
(as ‘lib’ is unspecified)
trying URL 'https://cloud.r-project.org/src/contrib/devtools_2.2.1.tar.gz'
Content type 'application/x-gzip' length 372273 bytes (363 KB)
==================================================
downloaded 363 KB

* installing *source* package ‘devtools’ ...
** package ‘devtools’ successfully unpacked and MD5 sums checked
** using staged installation
** R
** inst
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
*** copying figures
** building package indices
** installing vignettes
** testing if installed package can be loaded from temporary location
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (devtools)

The downloaded source packages are in
    ‘/tmp/RtmpfWDaCQ/downloaded_packages’

library(devtools)

Loading required package: usethis

install_github("WSpiller/RadialMR")

Skipping install of 'RadialMR' from a github remote, the SHA1 (18a7a6e9) has not changed since last install.
  Use `force = TRUE` to force installation

plot(by/bx, bx/byse, xlim=c(min((by-2*byse)/bx), max((bx+2*byse)/bx)), ylim=c(0, max(bx/byse)))
for (j in 1:length(bx)) {
lines(c((by[j]-1.96*byse[j])/bx[j], (by[j]+1.96*byse[j])/bx[j]),
c(bx[j]/byse[j], bx[j]/byse[j]))
}

abline(h=0, lwd=1); abline(v=0, lwd=1)

library(RadialMR)
r_input <- format_radial(bx, by, bxse, byse)

Missing SNP IDs; Generating placeholders

funnel_radial(egger_radial(r_input, alpha = .05))

Want to understand how all the pieces fit together? See the R for Data Science book: http://r4ds.had.co.nz/
Weights not specified: Adopting modified second-order weights

Radial MR-Egger

               Estimate  Std. Error    t value  Pr(>|t|)
(Intercept) -0.01214285 0.014905603 -0.8146503 0.4153372
Wj           0.44241552 0.005847025 75.6650678 0.0000000

Residual standard error: 0.728 on 2999 degrees of freedom

F-statistic: 5725.2 on 1 and 2999 DF, p-value: <2e-16
Q-Statistic for heterogeneity: 1588.522 on 2999 DF , p-value: 1

 Outliers detected